Large-Deviation Functions for Nonlinear Functionals of a 
Gaussian Stationary Markov Process 
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We introduce a general method, based on a mapping onto quantum mechanics, for investigating the 
large-T limit of the distribution P(r,T) of the nonlinear functional r[V] = (1/T) J* dT' V [X(T% 
where V(X) is an arbitrary function of the stationary Gaussian Markov process X(T). For T — > oo 
at fixed r we obtain P(r,T) ~ exp[— 9(r)T], where 8(r) is a large-deviation function. We present 
explicit results for a number of special cases including V(X) = XH(X) (where H(X) is the Heavisidc 
function), which is related to the cooling and the heating degree days relevant to weather derivatives. 
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I. INTRODUCTION 

The "persistence" of a continuous stochastic process 
has been the subject of considerable recent interest 
among both theoreticians and experimentalists in the 
field of nonequilibrium processes. Persistence is the prob- 
ability P(t) that a stochastic variable x(t) of zero mean 
does not change sign up to time t. Systems studied in- 
clude reaction-diffusion processes, phase-ordering kinet- 
ics, fluctuating interfaces and simple diffusion from ran- 
dom initial conditions [1]. Experimental measurements 
have been carried out on breath figures [2], liquid crys- 
tals [3], soap froths [4], and diffusion of Xe gas in one 
dimension [5]. In "coarsening" systems like these, which 
do not possess a definite length or time scale, the persis- 
tence has a power- law decay, P(t) ~ t~ e , at late times, 
and the persistence exponent 6 is in general non-trivial. 
In such systems, the normalized two-time correlation 
function, C(h,t 2 ) = (x(h)x(t 2 )) /[(x 2 (h)) (z 2 (i 2 ))] 1/2 
has the "scaling" form, C(t\,t 2 ) = f(ti/t 2 ), depending 
only on the ratio of the two times. In such systems, 
a simplification is achieved by introducing the logarith- 
mic time-scale T = hit, and the normalized variable 
X(T) = x(t)/ (x 2 (t)) 1 ! 2 , since the correlation function 
of X depends only on the time difference T 2 — T\, i.e. 
the process X(T) is stationary. Thus one is led to con- 
sider stationary stochastic processes. These processes 
are, of course, also of interest in their own right. In the 
new time variable, the persistence decays exponentially, 
P~ exp(-0T). 

In this paper we consider the integrated quantity 



dT'V[X(T')] 



(1) 



where V(X) is an arbitrary function of the stochastic 
variable X(T). Hence r is a functional of V. For the 
special case V(X) = H(X), where H(X) is the Heavi- 
side step function, r is just the fraction of time for which 



X(T') > in the time interval O < V < T. In this 
case, the probability distribution, P(r,T), of r for given 
T is just the "residence -time" distribution which, to- 
gether with the related "sign-time" distribution, where 
V(X) = sgn(X), has attracted a lot of attention earlier 
in the mathematics literature [6,7] and more recently in 
the statistical physics community [8] . Here we consider a 
general function V(X). We restrict our attention, how- 
ever, to the class of processes where X(T) is a station- 
ary Gaussian Markov process, for which some analytic 
progress can be made. A case of some interest, which 
we analyze in detail, is V(X) = XH(X). This case is 
relevant to weather derivatives as we now explain. 

It has been estimated that about 1 trillion dollars of 
the 7 trillion dollar US economy is weather related [9]. 
For example, weather conditions directly affect agricul- 
tural outputs and the demand for energy products and 
indirectly affect retail businesses [10]. Weather deriva- 
tives, first introduced in 1997, are financial instruments 
which allow hedging (by, for example, energy suppliers) 
against adverse weather conditions over a period of time. 
Here "adverse" might mean an unusually warm winter, 
when low demand for energy would affect the supplier's 
profits, as well as an unusually cold one when the supplier 
is unable to meet the demand. Temperature derivatives, 
the most common form of weather derivatives, are based 
on the concepts of "heating degree days" and "cooling 
degree days" , which are (rough) measures of the cumu- 
lative demand for heating and cooling respectively. 

Let X(T) be the temperature at time T in a given city. 
On a given day, n, the mean, X n , of the highest and 
lowest temperatures is recorded. The number of cooling 
degree days (CDD), over a period of N days, is given by 
CDD = Yl n =i max(X n — X Ql 0), where X Q is a reference, 
or baseline, temperature, while the number of heating de- 
gree days (HDD) is HDD = £^ =1 max(A - X n , 0). In 
the present paper we will, for simplicity, use an integral 
over continuous time rather than a sum, so that the cool- 
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ing degree days are given by CDD = J T dT' [X(T') — 
X ]H[X(T') - X ], where H(X) is the Heaviside step 
function. Thus CDD is the integrated temperature ex- 
cess (over the reference temperature) restricted to those 
periods where the temperature is above the reference 
level. It is a crude measure of the amount of cooling (air 
conditioning) required during the period T and also of the 
energy required to produce this cooling. Note that the 
power consumption of an air conditioner actually varies, 
for small temperature differences, as the square of the 
temperature difference between the room and ambient 
temperatures, so a better measure of the energy required 
would be J Q T dT' [X(T') - X } 2 H[X(T') - X ], instead 
of CDD. In the notation of this paper, CDD = rT, 
with V(X) = (X-X )H(X-X ) in Eq. (1). The num- 
ber of heating degree days (HDD) is similarly given by 
HDD = J T dT'[X a - X(T')} H[X a - X{T% which is a 
measure the amount of heating required in the period T 
and of the energy required to produce this heating. Esti- 
mating the likelihood of large deviations from the mean 
in CDD or HDD is clearly of interest to energy compa- 
nies. 

To make further progress, a realistic statistical model 
of temperature fluctuations is required. At this stage, 
however, any realistic model is too intractable to allow 
analytical progress. To illustrate the general approach we 
will instead employ a simple, though unrealistic, model in 
which the temperature, X(T), is taken to be a stationary 
Gaussian Markov process. We will discuss the limitations 
of this model, and possible improvements, in the Conclu- 
sion. As a further simplification, we take the reference 
temperature X Q to be the mean of X(T), though this 
simplification is not essential and can be relaxed. 

The outline of the paper is as follows. In section II we 
introduce the general approach to the problem of com- 
puting the distribution, Pr(r), of the quantity r defined 
by Eq. (1). Using a path- integral representation, the cal- 
culation of Pt{t) is mapped onto a problem in quantum 
mechanics, in which the function V(X) appears in the po- 
tential energy. In the limit of large T, only the quantum 
ground state energy is required. The final result takes 
the form Pt(t) ~ exp[— 0(r)T], where the function 9(r) 
is a large-deviation function. For the case of the sign- 
time distribution, corresponding to V(X) = sgn(X), r 
lies in the range — 1 < r < 1, and it is clear that 0(1) 
is just the usual persistence exponent, 9, since r = 1 re- 
quires X(T') > for < V < T. In section III, the 
method is illustrated on a number of special cases, of 
which V(X) = aX and V(X) = (3X 2 /2 are studied first 
as exactly-soluble tutorial illustrations before turning to 
the CDD problem (in the simple form outlined above) 
and finally the sign-time distribution. The last two ex- 
amples can be solved analytically in various regimes, and 
numerically elsewhere. Section IV concludes with a dis- 
cussion and summary of the results. 



II. GENERAL APPROACH 

Consider the general stationary Gaussian Markov pro- 
cess (Ornstein-Uhlenbeck process) dx/dt = —fix + 
where £(t) is Gaussian white noise with zero mean, and 
correlator (£(*)£(*')) = 2D5(t - t'). After the change of 
variables t = T//J,, x = (2D//j,) 1 ^ 2 X, the equation takes 
the form 



§ = -* + „(D. 

where (r](T)) = and 

(V(T) V (T')) = S(T T') 



(2) 



(3) 



The probability distribution of X(T') for < V < T 
is given by 



P[X] = Nexp 



dT' 



X 2 (T') + X 2 (T') } , (4) 



where X(T') = dX/dT', and N is a normalization con- 
stant. Our goal is to calculate the probability distribu- 
tion of r, defined by equation (1). In practice it is con- 
venient to look at the distribution P u (it) of the quantity 
u = rT. Its Laplace transform is 



P u (s) = (exp(-rsT)) = Z(s)/Z(0) , 
where Z(s) is given by the path integral 



(5) 



Z(s) = I DX(T) exp{-- J dT'(X 2 + X 2 

+2sV[X])}. (6) 

We are interested in the limit T — > oo. It is convenient 
to impose periodic boundary conditions, X(T) = X(0), 
since this restriction will not change the results in the 
large-T limit. Furthermore, the exponential in (4) should 
strictly contain the combination (X + X) 2 instead of 
(X 2 + X 2 ) . The missing term, 2XX, is however a perfect 
derivative, whose integral vanishes for periodic bound- 
ary conditions. Finally, with these boundary conditions 
the function Z(s) is the imaginary-time Feynman path- 
integral that gives the partition function of a quantum 
particle with Hamiltonian H = p 2 /2 + X 2 /2 + sV(X) at 
inverse temperature T, p being the canonical momentum 
conjugate to X. For T — > oo the ground state dominates 
to give, in this limit, 



(exp(-rsr)) = cxp{-T[£(,s) - £7(0)]} , 



(7) 



where E(s) is the ground-state energy for the Schrodingcr 
equation 



1 d 2 i\) 

Ydx 2 



U(X)i/> = E(s)i> , 



(8) 
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with potential 



U(X) = X 2 /2 + sV(X) . 



(9) 



For s = the problem reduces to a simple harmonic os- 
cillator, and E(0) = 1/2. 

The stochastic process x(t) studied above corresponds 
to the position of a Brownian particle in an external po- 
tential [ix 2 /2. For the case of a pure Brownian motion 
(fi = 0), Kac derived a formalism [6] to study the dis- 
tributions of arbitrary functional V[x] which also used 
a mapping to the Schrodingcr equation. Note, however, 
that the method presented above for the \i ^ case dif- 
fers in details from the original Kac formalism. 

To illustrate the method we discuss two simple exam- 
ples, before turning to some nontrivial cases, including 
the CDD problem. 



III. SPECIAL CASES 



A. V(X) = aX 



For this case we have 



rT = a [ dT'X(T') . 
Jo 



(10) 



This case is actually trivial since u, being a sum of zero- 
mean Gaussian variables, is itself a zero-mean Gaussian 
variable. All we require, therefore, is the variance, given 

by 



(u 2 )=a 2 dT, dT 2 <A-(Ti)A-(T 2 )) . (11) 
Jo Jo 

For the Ornstein-Uhlenbeck process (2), with noise 
correlator (3), one easily finds (X(Ti)X(T 2 )) = 
(l/2)exp(-|Ti - T 2 |). Inserting this result in (11), and 
extracting the leading large-T behavior, gives (u 2 ) —* 
a 2 T, and therefore (r 2 ) — ► a 2 /T. Hence the asymptotic 
distribution of r is given by (neglecting prefactors) 



P T (r) ~ exp(-Tr 2 /2a 2 ) , 



oo 



(12) 



We now show how the general machinery we have set 
up in section II recovers this result. The potential U(X) 
in the Schrodinger equation (8) takes the form 

U{x) = X 2 /2 + saX = (X + sa) 2 /2 - s 2 a 2 /2 . (13) 

This is just a harmonic oscillator with a shifted origin, 
so E(s) = 1/2 - s 2 a 2 /2 and, using (7), 



(exp(-rsT)) = exp(Ts 2 a 2 /2) 



(14) 



for large T. To recover Pr(r) we can invert the Laplace 
transform as follows. Neglecting pre-exponential factors, 



r ioc Hi 

P T (r) =T exp[T(rs + s 2 a 2 /2)\ . (15) 

J-ioc 2? ™ 

This integral can, of course, be evaluated exactly. Here, 
however, we use the method of steepest descents, which 
is valid for large T and can be readily generalized to the 
other, less trivial, cases that we will discuss. Writing 
the integrand in the form exp[Tg(s)], we have g(s) = 
rs + s 2 a 2 /2. The integral is dominated by the saddle 
point at s = —r/a 2 , where g = —r 2 /2a 2 . The integra- 
tion contour is deformed to pass over the saddle point, 
which lies on the real s axis. The saddle point is thus 
a minimum of g(s) with respect to variations of s along 
the real axis. The final result, ignoring non-exponential 
prefactors, is identical to (12). 
The structure of (12) is 



with 



P T (r) ~ exp[-0(r)T] , 



9(r) = r 2 /2a 2 . 



(16) 



(17) 



We can easily generalize this method to arbitrary 
V(X). The path- integral approach gives the general re- 
sult (neglecting prefactors) 



Mr) 



ds exp[Tg(s)] 



(18) 



where 



g(s) = rs + E(Q) - E{s) . (19) 
Using the steepest-descent method for the integral gives 



9{r) = max 



E(s) ~\-rs 



(20) 



where we have inserted E(0) = 1/2. The next example 
is a simple application of this idea. 



B. V(X) = (3X 2 /2 

This case illustrates the power of the method. The 
potential energy is now 



U(X) = (l + f3s)X 2 /2 



(21) 



so we have a harmonic oscillator again but with a mod- 
ified frequency, u = (1 + /3s) 1 / 2 , giving E(s) = (1 + 
13s) 1 ' 2 12. Thus 



9{r) = max 



i(l+/fo)Va_l_„ 




(22) 
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Note that 9(r) now has its minimum at r — /3/4, which 
is just the mean value of r [noting that (X 2 ) = 1/2 fol- 
lows from (2) and (3)], while large (r — > oo) and small 
(r — ► 0) values of r are strongly suppressed. An expan- 
sion of 8 near its minimum value gives a Gaussian dis- 
tribution P T (r) ~ cxp[-4T(r - /3/4) 2 //3 2 ]. This agrees 
with the result expected from the central limit theorem, 
i.e. Pr{r) ~ exp[— (r — /i) 2 /2er 2 ], with mean fi = /?/4 
and variance a 2 = /3 2 /8T. For T — > oo, the distribution 
becomes very narrow such that, at fixed r, the central 
limit theorem fails to give the correct asymptotics. The 
form (16) thus gives the behavior in the extreme tails of 
the distribution at large T. The function 8(r) is a "large- 
deviation function" which controls the distribution of r 
for large T. 

We note that this special case with V(X) = /3X 2 /2 
was also studied recently by Farago [11] in the context 
of power fluctuations in the Langcvin equation (2) by 
a somewhat different method. The probability density 
function of the dissipated power in reference [11] is pre- 
cisely the distribution Pr(r) studied here and the corre- 
sponding large-deviation function 6(r) has the same ex- 
pression as in Eq. (22). However our derivation seems 
simpler and easily generalizable to other forms of V(X) 
as we show below. 

C. V(X) =XH{X): "Cooling Degree Days" 

In this case the quantity rT gives the integrated value 
of X over the interval T, restricted to those values where 
X > 0. If X is the excess temperature over some base- 
line value where cooling becomes necessary, rT is a crude 
measure of the total energy consumption required to pro- 
vide the cooling [V(X) = X 2 H(X) would be a better 
measure, as discussed in the Introduction]. In the present 
simple model we have taken the mean temperature equal 
to the reference temperature (with X(T) being the de- 
viation from the mean), though this is not an essential 
restriction. The Schrodinger equation (8) takes the form 

i\>" - X 2 ^ - 2sXtP + 2E(s)i> = 0, X > , (23) 
ip" - X 2 ip + 2E{s)4> = 0, X < , (24) 

where ip" = d 2 ip/dX 2 . The required solutions can be ex- 
pressed in terms of parabolic cylinder functions, D p (x), 
using the standard solutions of the parabolic cylinder 
equation y" — (x 2 /4+ a) y = [12]. Selecting the solutions 
that satisfy the physical boundary condition ^>(±oo) = 
gives 

V-PO =*1>+{X) = AB p+ (V2{X + sj) , A>0, (25) 
= ip- (X) = B D p _ (-V2Jf) , X < , (26) 

where A, B are constants and 



P+ = E(s) - \ + *- , (27) 
P- = E(s) - \ . (28) 

The ratio A/B and the energy eigenvalues E(s) are ob- 
tained from matching the wavefunction and its deriva- 
tive at X = 0, i.e. we require ip+(0) = V'-(O) an d 
ip+'(0) — V'-'(O). This yields the eigenvalue equation 

IW(V2*) = D P -'(0) (2Q) 
V p+ (V2s) D p _(0) • 

The determination of E(s) from these equations is not 
possible analytically for general s. However, in the 
regimes s — > oo, and s — > — oo, which determine the 
small-r and large-r behavior respectively of 8(r), analyt- 
ical progress is possible. We consider these in turn. 



1 . The limit r — > oo 

As we shall see, in this limit it is sufficient to compute 
E(s) in the limit s — > — oo. This can either be done di- 
rectly from (29), or using the following (simpler) physical 
argument. 

Recall that, for V(X) = XH(X), the potential U(X) 
in the Schrodinger equation is 

U(X) = X 2 /2 + sX , X > , 

= X 2 /2 , X < . (30) 

For s — > — oo, the potential has a deep minimum, of 
depth — s 2 /2, located at X = — s. The wavefunction is 
exponentially small at X — 0, ip(0) ~ exp(— s 2 /2), and 
the change in the potential in the regime X < has 
an exponentially small effect on the ground state energy. 
Thus 

E(s) = -s 2 /2 + \/2 + 0{e- s2 ), s -> -oo . (31) 

The same result can be obtained, after some algebra, 
directly from (29). Inserting the result in (20) gives 
9(r) = max s (— s 2 /2 — rs) = r 2 /2. Since the maximum 
occurs at s = —r, the calculation is self-consistent for 
r — > oo. Thus we obtain 

6{r) r 2 /2 , r -» oo . (32) 

This large-r result has the same form as equation (17), 
which gives the general-r result for potential V(X) = aX 
(with a = 1 here). This correspondence is not so sur- 
prising: For r — > oo, the dominant processes X(T'), for 
< T' < T, will correspond to large positive X, so the 
step function H(X) in V(X) plays no role in this limit. 
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2. The limit r -> 

We shall see that this limit corresponds to s — ► oo. 
Again, one can use physical arguments as a short cut. 
For s -> oo, one has C/(X) = X 2 /2 for X < 0, with 
essentially a hard wall at X = 0. This gives, to leading 
order, E(s) — ► 3/2. We need, however, the leading cor- 
rection to this result. Since the wave function does not 
penetrate far into the wall, we can neglect the X 2 term 
in the potential for X > 0, i.e. write U(X) = sX. The 
Schrodinger equation then simplifies to 

il>" -2sXi; + 2E(s)ij = , X>0, (33) 
- X 2 V + 2£(s)V = , X < . (34) 

The wavefunction for X < is again a parabolic cylinder 
function, while for X > it can be expressed as an Airy 
function: 

V>(X) = AAi((2s) 1/,3 A) , A>0, (35) 
= BD p _(-V2A) , X<0, (36) 

where A, B are constants, and p_ = E(s) — 1/2 as be- 
fore. Matching the wavefunction and its first derivative 
at X — gives the eigenvalue equation 



(2g) 1 / 3 Ai , (0) 
\/2Ai(0) 



D P .'(0)_ ^ r (^ 

"d p _(o) r(-^) 



(37) 



where r(ar) is the Gamma function. In the limit s —> oo, 
the left-hand side of (37) tends to infinity, so the right- 
hand side (RHS) must also diverge in this limit. The 
ground state corresponds to the first divergence, where 
p_ — » 1. Therefore we write p_ = 1 — e in (37), 
and seek the leading behavior as e — ► 0. This gives 
RHS -> -(2/7r) 1 / 2 /e- Inserting this result in (37) gives, 
to leading order, 



= -2 1 / 6 



2 Ai(0)_ 1/3 



7T Ai'(0) 
2~2 1 / 6 r(l/3) 
7T 3V3 f(273) 



S -V3 = as -l/3 



(38) 



the last equation defining the constant a. Finally we have 
E(s) =p- + 1/2 = 3/2 - e. Inserting this in (20) gives, 
for r — > 0, 



6{r) = max(l - as" 1/3 - rs) = 1 - fer 1/4 



(39) 



where 



(?) 



3/4 



4^2 / r(l/3) \ 3/4 

3^ 3 / 8 vr(2/3)J 



2.04763... (40) 



Note that the value of s at which the maximum occurs 
in (39) is s = (a/3r) 3 / 4 , justifying our use of a large-s 
analysis of E(s) for the limit r — > 0. 



That 0(r) — » 1 for r — > is intuitively clear, since 
r = requires X(T) < for all T. This reduces to the 
usual persistence probability of the Markov process (2), 
for which 9 = 1. 



3. 9(r) for r near (r) 



Equations (32) and (39) give analytical results for 6(r) 
in the limits of large and small r respectively. For general 
r, 9 has to be computed numerically. There is, however, 
one other regime where analytical progress is possible, 
namely for r close to its mean value, where the central 
limit theorem (CLT) applies. 

The mean value is given by 



(r) = (V(X)) = (XH(X)) 
= I«X) + <|X|». 



(41) 



The Gaussian distribution for X gives immediately 
(X) = and (\X\) = ^/2/^(X 2 ) 1/2 = l/v^- Thus 



(r) 



2^ ' 



(42) 



In a similar way, the variance of r can be calculated, 
by exploiting the Gaussian properties of the process 
X(T). A tedious but straightforward calculation gives, 
for T -» oo, 



a 2 ee (O - (r) 2 = 



1 



2ttT 



(7T + hl2-2) 



(43) 



The central limit theorem then gives the behavior of 
Pr(r) for r near r c as PtM = (27rer 2 )~ 1//2 cxp[— (r — 
(r)) 2 /2er 2 ]. Inserting the explicit expressions for (r) and 
a 2 gives Pr(r) ~ exp[— 9(r)T], with 



(9(r) 



7T + In 2 - 2 



1 



2V^ 



(44) 



correct to leading nontrivial order in (r — (r) ) . 

The full result for 0(r) can be obtained by numerically 
solving (29) for the ground state energy E(s) for each 
value of r, then using (20) to find the corresponding 9{r). 
The result is displayed in Figure 1, with the asymptotic 
forms for r — > oo and r — > indicated. Note the very 
sharp rise to the value unity as r — > 0, as indicated by 
Eq. (39). 
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FIG. 1. The function 8(r) for the CDD problem, showing 
the asymptotic behavior for small and large r. 

In terms of the CDD problem, the behavior below the 
minimum (i.e. for r < (r)) determines the probability 
of an unusually small CDD ("cool summer"), while the 
region above the minimum corresponds to an unusually 
large CDD ("hot summer"). The fact that the function 
6(r) initially increases less rapidly to the right of the min- 
imum than to the left indicates that (within this very 
simple model) summers with a slightly larger than aver- 
age CDD are more probable than those with a slightly 
smaller than average CDD. This asymmetry is a conse- 
quence of the nonlinear relation between CDD and the 
temperature fluctuations, which are symmetric about the 
mean in our model. It should be stressed that the inte- 
gration period T has been taken to be large, to justify the 
steepest descent calculation. In practice this means that 
T (the length of a summer, say) must be large compared 
to the correlation time of the temperature (a few days, 
perhaps), which is not totally unreasonable. 



D. V(X) = sgn(X): The "Sign-Time Distribution" 

The structure of the calculation is similar to that of 
the preceding subsection. The Schrodinger equation is 



ip" - X 2 i> - 2[s - E(s)]il! = 0, X>0, 
i;" -X 2 tP + 2[s + E{s)]iP = 0, X<0. 

The solutions are parabolic cylinder functions, 

1>(X) =i>+(X) = AB P+ (V2X) , X>0, 
= ip-(X) = BD p _(-V2X) , X<0 

where now 

P± = E(s) T s — 1/2 . 



(45) 
(46) 



(47) 
(48) 



(49) 



Matching the wavefunction and its derivative at X = 
gives the eigenvalue equation 



D p+ '(0) = D p _'(0) 
D p+ (0) D p _(0) 



(50) 



which, using standard identities relating D p (0) and 
D/ (0) to gamma functions [12], reduces to 



r(-s±) ■ r(-^) 



(51) 



Although this equation cannot be solved analytically for 
general s, the limits s — > and s — > oo are tractable. 
Note that E(s) = E(—s) by symmetry, so there is no 
need to consider s — > — oo separately. 

The analysis starts from the potential well, U(X) — 
X 2 /2 + ssgn(X). For small s, the ground state energy 
E(s) is perturbatively close to, but slightly smaller than, 
£7(0) = 1/2. Therefore we write 



E(s) = 1/2 -6(a) , 



(52) 



where we anticipate S(s) — 0(s 2 ) from the symmetry 
E(s) — E(—s). Inserting this form in (49), (51) becomes 



(53) 



Expanding to fourth order in s and second order in 5 
gives 



S = s 2 In 2 -cs 4 + 0(s 6 ) , 



(54) 



where 



I 7T I 

c= -(In 2) 3 + —In 2 - -((3) = 0.269865.. . , (55) 

and C( n ) is t ne Riemann zeta function. Inserting this 
result in (20) gives 

6(r) = max[-rs - s 2 In 2 + cs 4 + 0(s 6 )] 



4 In 2 



+o(7 

V21n2/ V 



(56) 
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The maximum occurs at s = — r/21n2 + 0(r 3 ), so our 
study at small s is self-consistent at small r. 

For s — > oo, on the other hand, the potential develops 
a hard wall at the origin, and has a depth of — s next to 
the wall. Therefore we write E(s) = — s + 3/2 — e, with 
e small. Putting this form in (49), (51) becomes 

I> + e/2] r[e/2] 
r[ S + (e-l)/2] r[(e-l)/2] ' 1 j 

Taking the limits s> 1, and e<l readily leads to 

e = (tts)- 1 / 2 (58) 

to leading order for large s. Using this in (20) gives 

0(r) = max[l - (1 + r)s - l/(vrs) 1/2 ] 

S 

= l-- 2 (~d+r)) (59) 

The maximum occurs at s = [2y/n(l + r)}- 2 ^ 3 }, which 
tends to infinity as r — > —1 so the calculation is self- 
consistent in this limit. The symmetry of the problem 
under r — > — r leads to the more general result 

*(r) = 1 - | (^(1 - M)) ' + • • • , r - ±1 . (60) 

The function 9(f) is plotted in Figure 2, with only the 
region r > shown explicitly. The limiting behavior for 
small r and r close to unity is also shown. 



1 
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FIG. 2. The function 9(r) for the sign-time distribution, 
showing the limiting behavior for r — > and r — *• 1. 



IV. CONCLUSION 

In this paper we have presented a general method for 
computing the asymptotic behavior of the the probability 
distribution of the quantity r = (1/T) dT' V[X(T% 
where V(X) is an arbitrary function and X(T) is an 
Ornstein-Uhlenbeck stochastic process representing the 
motion of a Brownian particle in the presence of a stable 
parabolic potential /iX 2 /2. The main new result is that 
for fi > 0, the distribution of r, for large window size 
T has the form Pt(i") ~ exp[— 9(r)T], where 6(r) is the 
large-deviation function. The calculation proceeds via a 
mapping onto a quantum mechanical problem described 
by the Schrodinger equation (8) for a particle moving in 
the potential (9), where s is the Laplace variable conju- 
gate to rT. The inverse Laplace transform can be per- 
formed using steepest descents in the limit T — > oo. 

Two non-trivial applications have been presented. The 
first is a calculation of the large-deviation function 9(r) 
for the cooling degree days problem, based on a sim- 
ple model of temperature fluctuations. The model used 
assumes that the temperature fluctuation X(T) is de- 
scribed by the Ornstein-Uhlenbeck process (2), i.e. that 
X(T) is a stationary Gaussian Markov process with time- 
independent noise strength. This may not be a realis- 
tic model for several reasons: (i) A simple Markov pro- 
cess is not thought to be an optimized model of tem- 
perature fluctuations, which tend to exhibit stronger au- 
tocorrelations than a Markov process. A more plausi- 
ble statistical model, according to [10], writes the fluc- 
tuation, X n , of the mean (average of daily high and 
low) temperature on day n as the linear combination 
X n = Y^m=\ w n.-raX m + r] n , where the "memory ker- 
nel" w n - m is a decreasing function of n — m and rj n is 
uncorrelated Gaussian noise. The Markov case corre- 
sponds to k = 1. Here X n is the difference between the 
measured mean temperature on day n and its expected 
value. The latter should contain a 365 day seasonal vari- 
ation (roughly sinusoidal) (ii) The noise strength should 
also contain a 365 day seasonal variation: the variance 
of the temperature fluctuations can be different at dif- 
ferent times of the year, (iii) The reference temperature 
for the calculation of CDD and HDD should in general 
be different from the expected temperature. We hope to 
incorporate some of these features in future studies. 

The second non-trivial application is to the "sign-time" 
distribution. In the context of "sign-time" , we point out 
that the asymptotic form of the "sign-time" distribu- 
tion has a markedly different behavior in the Ornstein- 
Uhlenbeck process (where a Brownian particle moves in 
a stable parabolic potential fiX 2 /2) studied here com- 
pared to the ordinary Brownian motion (fj, = 0). In 
the later case, the "sign-time" distribution Pt(t) is in- 
dependent of the window size T for all T and is given by 
Pr(r) — \ jn\J\ ~ r 2 [7]. In the former case (ji > 0), on 
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the other hand, the "sign-time" distribution depends ex- 
plicitly on the window size T and Pt{t) ~ exp[— 9{r)T] 
for large T where the large-deviation function 6{r) has 
been computed exactly in this paper. 

We further note that for this "sign-time" problem, the 
function 6(r) can also be obtained using the "independent 
interval approximation" (IIA) [13,14], which exploits the 
fact that the intervals between zero crossings arc sta- 
tistically independent for a Markov process. In fact, 
for renewal type processes where the successive intervals 
are statistically independent, the "sign-time" distribu- 
tion has been computed by Godreche and Luck [15] using 
the interval size distribution as an input. Their result can 
be simply extended to calculate the "sign-time' distribu- 
tion for the Ornstein-Uhlenbeck process. The IIA also 
has the virtue that it can be used to obtain approximate 
results for non-Markov processes. Persistence exponents, 
for example, are often given rather accurately by the IIA 
[1]. However, it is not straightforward to adapt this IIA 
method to general nonlinear functions V(X), whereas the 
path-integral approach and mapping onto quantum me- 
chanics adopted here is readily applicable to any V(X). 

Here we have only considered Gaussian Markov pro- 
cesses mainly because they are simple and amenable to 
analytical calculations. Recently the calculations of the 
asymptotic distributions for the "sign-time" and other re- 
lated quantities such as "local-time" have been extended 
to non-Gaussian Markov processes where a Brownian 
particle moves in an arbitrary stable or unstable potential 
and moreover exact results have been obtained [16] when 
the underlying potential is random as in the Sinai model. 
The extension of these methods and results presented 
here to non-Markov processes, however, still remains as 
one of the outstanding challenges for the future. 
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